Characterization of diffractive bifocal intraocular lenses

Multifocal intraocular lenses incorporate a variety of design considerations, including dimensioning of the base monofocal shape and the diffraction grating. While studying three different lens models, we present a practical approach for mathematical modelling and evaluation of these geometries. Contrary to typical lens measurement methods, non-contact measurements were performed on the Alcon SN6AD1, HumanOptics MS 612 DAY and the AMO ZMA00 lenses using a confocal microscope. Subsequent data processing includes centering, tilting correction, filtering and an algorithmic decomposition into a conic and polynomial part and the diffraction grating. Lastly, evaluation of fitting parameters and grating shape is done to allow for inferences about further optical properties. Results and analysis show the confocal microscope to be a suitable imaging method for lens measurements. The processing of this data enables the reconstruction of the annular diffraction grating over the complete lens diameter. Apodization, near addition and diffraction efficiency characteristics are found utilizing the grating shape. Additionally, near-optical axis curvature, asphericity and higher order polynomials are identified qualitatively from the reconstruction of the monofocal base form. Derived properties also include the lens optical base and addition power. By making use of the surface geometries, as well as the lens’ material and thickness, a full lens model can be created for further studies. In summary, our analytical approach enables the insight to various intraocular lens design decisions. Furthermore, this procedure is suitable for lens model creation for research and simulation.


Materials and methods
All processing is implemented in Python, while utilizing the well-known scientific libraries numpy, scipy and matplotlib. The code is hosted in the following repository: https:// github. com/ droch eam/ miol-reng-tools.
Lenses. Investigated lenses include three different models from three different manufacturers, two different materials and three different base and addition powers. While different design aspects would be better compared for models with equal powers, varying values and therefore highly different surface shapes enable us to test the characterization procedure more thoroughly. The examined lenses are as follows: Alcon SN6AD1. The examined AcrySof ® IQ ReSTOR ® SN6AD1 intraocular lens is specified with 13.0 dpt of optical power and +3.0 dpt of near addition. It is made of a highly refractive hydrophobic acrylic material with a refraction index of n = 1.55 . The overall shape is symmetrical biconvex and introduces −0.1 µm of spherical aberration. The apodized diffractive structure is located on the anterior aspheric surface 6,7 .
HumanOptics MS 612 DAY. Secondly comes the HumanOptics (formerly Dr. Schmidt) MicroSil ® 612 DAY intraocular lens with 25.5 dpt optical power and +3.5 dpt near addition. It is made of silicone material and is built as biconvex shape with aspheric surfaces, whereas the anterior side holds the diffractive profile 8,9 . Unfortunately, no data on the refraction index was available. Subsequently a value of 1.47 is assumed, which lies within the typical range for silicone.
AMO ZMA00. The third lens is an AMO Tecnis ® ZMA00 lens with 30.0 dpt optical power and a +4.0 dpt addition. It consists of hydrophobic soft acrylic material with a refractive index of n = 1.47 . The lens is built as a pupil independent, fully diffractive multifocal posterior surface and an aspheric surface on the anterior side. −0.27 µm of spherical aberration are introduced by the lens geometry 2,6,10 .

Measurement.
Commonly used methods for lens measurement include interferometry and profilometry setups 11 . Most of these are highly specialized or cost-ineffective. This work utilizes an available µsurf ® custom confocal microscope by NanoFocus AG, Oberhausen, Germany. In contrast to contact profilometry, the resolution is independent of the stylus size, moreover the lens surface can not be damaged by contact. Additionally, the two-dimensional data output allows for subsequent centering, noise removal and surface smoothing.
Measuring the different microscope focus distances on the edge of the optical relevant area of the lens makes it possible to calculate the tilt vector of the surface. The tilt values are then compensated with the help of a micrometer platform, whereas a residual tilt is removed using software processing shown in "Alignment and centering". The entire height data consists of multiple sub-images, which need to be stitched together to form one continuous height image. But for the stitching to work effectively, we need to enhance the surface's information content. This is done by introducing artificial surface impurities to the otherwise smooth lens topography. These consist of microscopic graphite particles with a height around 0.5µm, which will be visible as contamination in the measured data. Sufficient graphite contamination enhances the success rate of the shift detection severely.
The measurement consists of a row of 36 images, with each image being 320 µ m × 320 µ m in size. To improve shift detection, the overlap region between images was set to as much as 150µm. In regard to resolution, in the lateral dimension it is quantified to 0.625µm, while the axial resolution of each image is expected to be better than 0.1µm.
With the confocal measurement being a non-contact method, there are no limitations on the lens material. However, all surfaces need to be in a dried condition, since every remaining liquid film will be interpreted as part of the surface. Also, hydrophilic lenses are known to slightly change their shape while drying, but we currently don't have any empirical data on how severe the impact on measurements and results will be. An important requirement however is that the surface needs to be reflective for the wavelength range of the microscope illumination. Another aspect is the maximum measurable surface slope, which is limited by the numerical aperture of the objective. For higher slopes the light is reflected away from the beam path of the microscope. In our measurement setup the angle range is limited to ±26.6 • , which is sufficient for the examined surfaces.
Image shift detection. The stepper motor of the microscope is responsible for displacement between individual image acquisitions. However, with the motor being an open loop system, there is no precise feedback on its absolute position and consequently the lateral image shift values. The typical image position uncertainty lies in the order of a few micrometers, but adds up over a large amount of images. With the help of overlapping image regions and methods for shift detection, the shift vectors can be estimated.
Contrary to the built-in shift detection method of the microscope, we implemented a different technique that can be fine-tuned to our needs. This algorithm utilizes the Fourier transformation's shift theorem. Let the height functions in the overlapping region be with the assumption of the second image being the first one shifted by a vector � s = (s x , s y ) . x and y denote the lateral dimensions, while the axial dimension of the microscope is called h. Transforming to the Fourier domain yields: www.nature.com/scientificreports/ The exponential function Q is isolated using the ratio Inverse transformation then yields a Dirac pulse δ at the image shift position.
Finding the maximum in V obtains the position of the Dirac pulse and therefore the shift vector s . Parameters for the stated algorithm include the search region and a threshold for the detection of the peak. In practice, Eq.
(1) does not always hold true, because of noise and some information only being present in one of the pictures. Furthermore, the algorithm needs sufficient surface details to find a shift vector. In multiple lens measurements this approach identified shift values for more than 90% of all image transitions, the rest being explained by missing graphite particles in those regions.
Alignment and centering. The lens center and diameter are elementary for the further processing steps.
Within the processing script the user specifies the lens center and edge with the help of a surface derivative image, where slight changes, including the grating rings, are distinctly visible. Next, an optimizing algorithm determines the lens tilt by maximizing rotational symmetry around the lens center. Thereafter, the script corrects the tilt by applying these calculated values.
Profile decomposition. With a radially symmetric surface, a one-dimensional profile is sufficient for describing the whole topography. For that, the processing script converts the two-dimensional surface image into an average, one-dimensional profile reaching from lens center to the outer edge.
Next, the profile is decomposed into multiple parts for analytical modelling: The radially dependent lens profile z consists of a conic section part z c , a polynomial component z p and the diffraction grating z d , as specified by Eq. (5).
Fitting of base shape z c (r) + z p (r) and grating z d (r) are performed independently by assuming that the local changes due to z d don't influence the base form consisting of z c and z p . The z c component contains an offset z 0 and the well-known conic section formula.
The aim is to determine curvature ρ , conic constant k and the offset z 0 best matching the acquired data. For that reason Eq. (6) is rearranged to: Equation (7) can be represented in the form of a linear equation system b = Ax with: where r 1 . . . r i and z 1 . . . z i are the measured values. The overdetermined equation system is best solved using the QR eigenvalue algorithm, which outputs a least-squares fit for x at a low computational cost. From x the desired values z 0 , ρ and k are found. To minimize the influence of the polynomial component, which is also part of z(r), only data in the inner 75% of the lens radius is used, where z p plays a minor role.
Polynomial regression on the difference z(r) − z c (r) with even orders up to n = 10 determines the polynomial component z p . According to Eq. (5) the remaining diffractive part is then z d (r) = z(r) − z c (r) − z p (r) . In the last step the diffractive profile is fitted using sectionwise polynomial functions of fourth order.

Results
Base shape. Results of the decomposition are shown in Table 1 for the anterior side and in Table 2 for the posterior lens side. R and k denote the surface curvature and conic constant derived using the fitting in "Profile decomposition", where R is the inverse of the curvature parameter ρ . For the sake of simplicity, the height change h p due to a polynomial component is listed, instead of specifying individual coefficients. h denotes the total height. d o is the optical diameter, thus twice the highest radial distance on the surface. Each value is the mean of two independent measurements for every lens side. www.nature.com/scientificreports/ Table 3 shows the resulting properties of the lens with n being the refractive index known from manufacturer data, d e the size of the lens edge, d o being the difference in optical diameter of back and front side. The overall thickness d is the sum of h of both lens sides as well as the edge thickness d e . d e was determined with an additional confocal measurement of the lens edge. D denotes the base power of the lens and D add the near addition. A graphical overview of the geometrical quantities is illustrated in Fig. 1.
The optical base power results from the lens maker Eq. (9) with the aqueous humour index of n a = 1.336 15 . Further, the near addition power D add is the result of calculations described in "Near addition power".
Although being declared as symmetrical biconvex, the SN6AD1 shows slightly varying curvature radii R 1 and R 2 , the same being the case for the HumanOptics MS 612. The AMO ZMA00 has strongly varying radii, with a ratio of around 1:2.2 between anterior and posterior side. The Alcon SN6AD1 and AMO ZMA00 show a conic constant k < −1 on both lens sides, resulting in an outwards declining surface curvature and therefore a lower optical power outside the lens center. This is the expected behavior for lenses with negative spherical aberration, as is specified by the data sheets of the lenses. On the other hand, conic constants k 1 and k 2 varying in sign, like for the MS 612, are an indication for a lens design with no added spherical aberration. www.nature.com/scientificreports/ Lenses MS 612 and ZMA00 feature a higher order polynomial component h p , while the minor polynomial parts for the SN6AD1 front side are probably due to measurement uncertainties and processing artefacts. The thickness d of the HumanOptics and AMO lens is similar, whereas the Alcon model has a thin lens design with only d = 0.43mm. This is the result of the lower optical base power, which produces smaller heights h 1 and h 2 .
Another interesting aspect is the optical diameter difference d o in Table 3. The front of the SN6AD1 is smaller by roughly 90 µm and for the ZMA00 by 110µm, respectively. This design choice is motivated by minimizing stray light: An angled lens edge reflects impinging light away from the inner parts of the retina, a curved lens edge additionally distributes the light heterogeneously inside the eye 12 . Evidently, Alcon and AMO are known to incorporate such a design in similar lens models and families 13 .
The curves coincide with the expected behavior of a kinoform phase grating. However, the innermost zone of the AMO ZMA00 deviates by consisting of two linear segments. The SN6AD1 grating shows eight steps, the MS 612 nine steps and the ZMA00 lens 29 zone steps. While the diffractive profile follows through the whole surface of the ZMA00, for the two other lenses it only exists at the inner 3.6 mm diameter.   www.nature.com/scientificreports/ The Alcon SN6AD1 shows a distinct decreasing apodization, whereas the zone edges have nearly constant height for the AMO ZMA00. While a varying step height is visible on the HumanOptics MS 612, part of it could be a visualization artifact as discussed later in "Design parameter deviation".
In the lenses' mean profile curves there are still some remains of noise and surface contaminations like dust or graphite particles, especially noticeable in the HumanOptics MS 612 anterior profile. Aside from that, the profiles in Figs. 3 and 4 show a superimposed waviness with an amplitude of around 500 nm, which can be a result of stitching errors or polynomial fitting artifacts.
Near addition power. The optical power of a kinoform profile is directly dependent on the annular zone positions r i . Figure 5 illustrates such a kinoform grating with its focus, whereas the spherical wave fronts moving towards the focus are traced back to the grating to find grating intersections. One can see a dependency of r i on the zone number i = 1, 2, . . . , the design wavelength 0 , the focal distance f add as well as the offset z at the optical axis.
The offset at the optical axis z can also be described using a phase offset φ 0 : From the right-angled triangle between the origin, position r i and the focus follows: Solving for r 2 i yields Eq. (12).
With φ 0 ∈ [0, 2π) the ratio φ 0 2π is bound to [0, 1). Also, for a typical lens f add ≫ 0 and f add ≫ 0 (i − 1) 2 can be assumed. This simplifies Eq. (12) to: φ 0 is commonly set to 2π , resulting in: This relationship is utilized for the MS 612 and ZMA00 lenses. Alcon on the other hand sets φ 0 = π , giving us Eq. (15) for the SN6AD1 lens 4 . This varying grating offset is already visible in Fig. 2, where the innermost part at r = 0 µm starts at half the zone's step height, compared to the full height seen for the two other lenses. www.nature.com/scientificreports/ Solving Eqs. (14) or (15) for f add makes it possible to determine the power D add for each phase zone. Mean values for all zones on each lens are found in Table 3. A comparison of model and measured zone edge locations is displayed in Fig. 6. The lenses SN6AD1 and MS 612 show satisfying compliance with the expected behavior, while the radii of the ZMA00 show a consistent downward deviation.
Design parameter. The design parameter α i is defined as 14 : With n being the material's refractive index, n a the aqueous refractive index, 0 = 550 nm the design wavelength and h i the height of the i-th zone step. The design parameter directly influences the diffraction efficiency η j for the j-th order, as described by Eq. (17) 14 .
Design parameter values α 1 derived from the innermost zone are depicted in Table 4. Parameters α 0 for both SN6AD1 and ZMA00 are known from literature. However, α 0 is specified for the center of the lens, deviating from α 1 for apodized lenses as the SN6AD1. Nevertheless, the measured values are in proximity to the expected parameters, but show a downward deviation.

Discussion
Optical power deviation. Compared to manufacturer data in "Lenses", the optical base powers D in Table 3 show small differences. In terms of value this amounts to 0.25% relative error for the SN6AD1, −0.12% for the MS 612 and the ZMA00 having the highest error with 1.53% . All results were within the manufacturing tolerances according to EN ISO 11979-2: 2014 15 . The small metrological deviations are therefore understandable and would have no impact in everyday clinical practice. An overview of designed D, D add and measured quantities D m , D add,m , as well as their relative errors D rel , D add,rel can be found in Table 5. According to Eq. (9) deviations in optical power may arise from variations in quantities n, R or d. The curvature radii R are expected to be fairly accurate, since the total surface shape contributes to this quantity. Due to a tilted edge and measure-  www.nature.com/scientificreports/ ment uncertainties of both lens surfaces, we don't expect the total thickness d to be more exact than d = ±40 µm. However, the thickness change d would account for less than 0.01dpt following Eq. (9) in all investigated lenses. The refractive index n is only specified to two decimal places, the index difference (n − n a ) is therefore accurate to merely two significant places. This possibly explains most of the deviation at hand. Nevertheless, for being derived from measured data and simplified considerations, the results provide an excellent fit to manufacturer data. Relative errors of the near addition D add from Table 3 are 1.92% for the SN6AD1, 1.02% for the MS 612 and a value of −5.94% for the ZMA00. A possible explanation is the location of the diffractive part: While not lying in the principal plane of the lens, its optical power needs to deviate to produce the same effective power. One would expect a lower power value for an anterior diffractive side and a higher value for a posterior one. As it turns out, this exactly matches the deviation's direction for all three intraocular lenses. However, the discrepancy would be best resolved by simulation of these geometrical models in further research.
Design parameter deviation. Differences in the design parameter α in Table 4 are more severe. A change of α = 0.51 to α = 0.47 produces an absolute efficiency decrease �η 0 of −6.5% for the zeroth order. While this form of profile visualization is suitable for a rough estimate, it is not suited for an exact characterization. One possible cause for the deviation is the non-zero width of the zone edges, resulting in a finite edge slope. The grating edges show a typical width of 2-5 µ m in the radial dimension. This width is based on manufacturing artifacts, although it can be a design decision as well. It is known, that the zone edges result in optical dead zones, regions where rays are refracted a second time, producing stray light and leading to an efficiency loss. By introducing a so-called groove angle, a specifically sloped zone edge, the effect can be compensated 16 . This solution is not only found in Fresnel lenses, but Zeiss is known for incorporating it in the ZEISS AT LARA 829 MP intraocular lens 17 .
The step width leads to a decrease in the step height, since the transition between one zone to another starts at an earlier and ends at a later radial position, disrupting the initial form of the annular zones. Another factor is the converging and pre-focused light coming from the cornea. An angled ray "sees" the projection of the step height, therefore either the height or the slope angle of the steps needs to change to keep the projected height constant. In the latter case an increasing groove angle should be incorporated in the grating design, since the mean incident angle also increases radially. Furthermore, the apodized profile of the SN6AD1 leads to an additionally decreased step height, when measured anywhere except at the lens center.
Finally, some missing height can be explained by the filtering and fitting artifacts from the software processing.
Comparison to earlier work. Diffraction profile properties of the ZMB00, an MIOL similar to the ZMA00, were already examined by Loicq et al. 2 . But compared to results of our full field measurements only a minor profile section is presented there. Nevertheless, the step heights are showing a similar magnitude. Furthermore, the zone edge positions r i also show a downward deviation from the model values, although it is less pronounced. However, a lens with base power of 20dpt is measured, while our model has 30dpt. The lens base power has an effect on the lens thickness and therefore the position of the grating relative to the principal plane of the lens. Hence, previous considerations from "Optical power deviation" could explain the higher deviation for this lens. The surface shape of the monofocal version of the Tecnis ® , the ZCB00, was studied by Miret et al. 3 . Their work arrives at the conclusion of the anterior side being conic and the posterior side being spherical, both sides without higher order polynomials. These results differ from ours, where both surfaces are aspherical with a polynomial component. It is not obvious from our research, whether the deviation arises from an incorrect surface component decomposition or from the shape of the multifocal version actually being different.
In the case of the Alcon SN6AD1 it is known from Madrid-Costa et al. 18 , that the grating consists of eight steps, therefore nine zones, inside the inner 3.6 mm diameter of the lens. Additionally, the grating starts at around half the step height of the first zone. This is in accordance to the results of the work at hand.
Symmetrical, yet slightly different curvature radii of posterior and anterior surface, as well as the highly negative conic constant of the front, are known from manufacturer data of the similar monofocal lens model SN60WF 19 . Unfortunately, no direct comparison is viable, since there is no data provided for a model with the same optical power. Data quality. Regarding data quality a differentiation between overall surface shape and derived surface components must be made. For instance, the conic constant k is tainted with larger uncertainties, since the profile differs barely from a sphere in the investigated area. The fitting is therefore largely impacted by stitching errors or surface differences mistakenly being attributed to the polynomial component. The same is true for the polynomial part, which proportionally accounts for even a smaller height change. Attribution of surface behavior to an incorrect component does not affect the overall lens shape, since it includes all components altogether. Furthermore, although the superimposed waviness of about 500 nm in the diffractive part in Fig. 3 appears large, it is minor compared to the whole height change of h 1 = 470 µm of the surface shape.
Nevertheless, it has to be acknowledged that the decomposition technique from "Profile decomposition" comes with its difficulties. Besides trouble distinguishing surface behavior between the polynomial and conic component, the conic fitting region is a source of result variation. Increasing the fitting region further from 75% of the diameter would lead to slightly different surface components. Another issue is filtering and fitting, due to their low pass properties some abrupt changes in the grating are potentially smoothed out. Further improvements could be made with the usage of a microscope with a precise closed loop stepper motor, eliminating the need for and the errors from the shift detection algorithm. www.nature.com/scientificreports/ An alternative asphere fitting approach, including conic section and higher order polynomials, is known from Sun et al. 20 . Their solution is based on an initial linear least squares estimation and an additional non-linear leastsquares step. It is not clear, if the mentioned approach would be suited for a surface with an additional grating and grating center height offset like in our application.
Altogether, the data quality is sufficient to identify important properties, as could be seen in the sections before.

Outlook
The obtained mathematical model should provide an excellent starting point for simulation using established raytracing software like Zemax OpticStudio ® or OSLO ® from Lambda Research Group. Simulations of these lens models, especially in combination with an adequate eye model, will show how close the model really comes to the expected behavior of the product. The monofocal base shape of the lens can be modelled using the parameters from Tables 1 and 2 and an aspheric surface description in both programs. In contrast, our script exports the whole shape including the grating as a spacing and height dataset. While both programs support user defined surfaces, providing details on how to import this data is outside our expertise.
In further research the processing procedure could be extended to support analysis of trifocal or EDoF (Extended Depth of Field) IOL. While there would be no differences in measuring those lenses, the analysis and mathematical modelling require a different approach. This is especially true for surface shapes that have no rotational symmetry or can't be decomposed into a simple base and differential profile.
Finally, the presented approach could be performed on more up-to-date and clinically relevant multifocal intraocular lenses to allow for detailed insight into the current lens designs.